Investigation the effect of water addition on intermolecular interactions of fatty acids-based deep eutectic solvents by molecular dynamics simulations

In this work, we focused on the interaction between hydrogen bond acceptor (HBA) and hydrogen bond doner (HBD) in the binary mixtures. The results showed that Cl− anion plays a key role in the formation of DESs. Also, the structural stability of deep eutectic solvents based on fatty acids (FAs) and choline chloride (Ch+Cl−) at different ratios was investigated in water using molecular dynamics simulations. We observed that the interaction between the chloride anion and the hydroxyl group of the cation leads to the transition of HBA to the water-rich phase. These atomic sites have important rule in the stability of the eutectic mixtures based on FAs and Cl− anion. However, it seems that the binary mixtures with the mole percent at 30% of [Ch+Cl−] and 70% of FAs have more stability than other ratios.

www.nature.com/scientificreports/ and acceptors in hydrophilic DESs. However, since these solvents have low toxicity, a feature that makes them suitable for the pharmaceutical industry. Despite various applications of hydrophilic DESs (high miscibility), this feature limits the utility of hydrophilic DESs for polar solutions 8 . However, DES-based on hydrophilic compounds have unique benefit such as the desulfurization of liquid fuel and natural gas and synthesis of nanomaterials, and electrochemistry. In the design of hydrophilic solvents utilizes choline chloride salt and quaternary ammonium salts (QAS) as HBA. However, the length of the alkyl chain components of DESs could play a crucial role in the design of stable hydrophobic solvents. Hydrophobic deep eutectic solvents were introduced by van Osch and co-workers in 2015 8 . Hydrophobic DESs have a lot more efficiency in the removal of pollutants from aqueous. Accordingly, van Osch and co-workers showed that HDES-based quaternary ammonium salts and decanoic acid (DecA) have high extraction yields in the separation of volatile organic compounds 9 .
There is great research to investigate the stability of DESs in water by using MD simulation. A molecularlevel understanding DESs in the presence of water and further investigations are necessary in this direction. Further study by the authors of this work was done in order to investigate the relationship between the percentage composition of the two components of HBA and HBD and the stability of DES in adjacency water. Partial disruption of the structural properties of the DESs occurs when DES is soluble in water. The fatty acid and choline chloride-based DES integrity is lost to some extent in adjacency water. The addition of the water molecules to box simulation containing deep eutectic solvents generally causes the dynamics to become fast. In this work, through molecular dynamics simulations, the influence of water on the dynamic and structural properties of DESs based on fatty acids (Caprylic acid, Lauric acid) and choline chloride were investigated in the adjacent water. The binary mixtures at different molar ratios of Ch + Cl − and acid were simulated in the adjacent water and the pure state at 353 K. In order to check the stability of DESs in the adjacent water, the distribution of two-component of DESs around each other was investigated by calculating the combined distribution functions (CDFs), radial distribution functions (RDFs), and angular distribution functions (ADFs), and spatial distribution functions (SDF). Furthermore, the dynamics of DESs is characterized by computing the self-diffusion coefficients, the mean square displacements (MSDs) and velocity autocorrelation functions (VACFs) of the center of mass (COM) of species as a function of time t . We find that the binary mixtures consisting of DES1 (Caprylic acid: Choline chloride) and DES2 (Lauric acid: Choline chloride) have stability in adjacency water.

Methodology
Initially, the structures of eutectic solvent componentes and pollutant molecules were drawn in Visual Molecular Dynamics (VMD) software. Then structures were optimized at the MP2/6-31G* theory level using the Gaussian 09 software 10 . Force field parameters and partial charges of species were obtained by fitting results from quantumchemical calculations (More details about parameterization are given in the previous work) 11 . The results of this paper present that the shear viscosities obtained for deep eutectic solvents based on caprylic acid and thymol at a ratio of 1:1 using MD simulation are close to the experimental value, with a difference of less than 8.82% 11 . To preparation of initial simulation boxes, packmol-16.343.3 package were employed (The pachage is distributed as free software and can be downloaded from http:// www. ime. unica mp. br) 12 . In the beginning, the cubic box containing 700 choline chloride ion pair and 300 fatty acid was randomly provided. Then, the mixtures are composed of species with the composition of fatty acid/choline chloride being 500: 500 and 700:300. Periodic boundary conditions were imposed in all directions to mimic bulk properties. Molecular dynamics simulation for all the studied systems was performed with time step 1 fs using the NAMD-2.12 package (https:// www. ks. uiuc. edu/) 13 . The information of the simulated boxes and the different fatty acids and choline chloride salt with abbreviations are listed in Table 1. The structure of the compounds of the binary mixtures is shown in Fig. 1. First, the 5,000,000-step energy minimization was undertaken to remove bad van der Waals contacts of the binary mixtures. The binary mixtures were heated to a temperature of 353 K. All binary mixtures were equilibrated via a 50 ns in the NPT ensemble. The Langevin piston Nose-Hoover method and Langevin dynamics were applied to keep at constant temperature of 353 K and pressure, respectively. The equations of motion were solved by using the standard Verlet algorithm with a time step of 1 fs. The particle-mesh Ewald (PME) algorithm was used to calculate electrostatic interactions and a cutoff of 12 Å was set for Lennard-Jones interactions 14 . After the system reaches equilibrium, the dynamic and transport properties were investigated and the last 1 ns of the simulations was considered to combine with the water box. The two cubic boxes were combined together to get the periodic cell that we will refer to this as the binary systems in the adjacent water. The binary systems in the adjacent water were investigated under the same conditions of pure state.
Radial distribution and combined distribution functions. The radial distribution function, g(r), and coordination numbers (CN) give information concerning the structural properties of the binary mixtures. The RDF, g(r), is related to the coordination number, or the number of neighbors, by the following equation: www.nature.com/scientificreports/ where the bulk density denoted by ρ. g(r) is also called the pair distribution function and can describe the distribution of the species around any given specie in the system 15 . The site-site radial distribution functions (RDFs) have led researchers to examine in detail the structural properties. Therefore, special attention will be paid to site-site RDF analysis in this section. The different atomic sites of the species for drawing RDFs are the oxygen and hydrogen atoms of carboxyl (-COOH) of FAs and hydroxyl (OH-) functional groups of the Ch + cations. This analysis can be useful for H-bond analysis, which will be examined in detail. It is the HA and OA atoms of FAs and the O and HA atoms of Ch + cation. RDFs between HA atom of FAs and Cl − anion in pure binary mixtures, as shown by g(r) Cl − −FAs , represnt a significantly sharp peak compared to other atoms of FAs at around 2.0 Å (see Fig. 2a). In addition, the peaks of all the RDF between the Cl − and different atoms of HBD molecules appeared at long distances. Higher occurrence probabilities of H-bonds can be observed at 2.0 Å for the Cl − −HA FAs interaction. Angular/distance probability region at 105°−120°/2.0−2.5 Å is a confirmation of the RDF results (see Fig. 3a). The structural correlations of HBA-HBD and the orientation of the individual solvent molecules were analyzed by combined distribution functions (CDFs). Angular Distribution Functions (ADF) were calculated using TRAVIS-200504 (This pachage can be downloaded from http:// www. travis-analy zer. de/) 16 . ADFs of molecular dipole vectors are ontained from as: The position vectors can be used for specifying the angle of interest 17 . The preferred orientation was investigated by combining the RDF between the HA atom of [Ch + ] and the O2 atom of FAs and an ADF between the C1 FAs −O2 FAs −HA Ch+ angle. Also, the angle between the C1-O1 vector in FAs and the intermolecular Cl − -O1 FAs vector between anion and FAs can be defined for anion orientation around FAs. There is the maximum occurrence probability region at an angle range of 105−120° in the less distance ∼3 Å for anion and FAs interaction in the binary mixtures (see Fig. 3a). The maximum occurrence probability region was observed for Ch + cation and FAs at around 150°-180°/2 Å (see Fig. 3b). CDF analysis in the different ratios confirms that the  Fig. S1a-c). The g(r) Cl − −FAs and RDFs between choline cations and FAs represent the sharp peak at 2 Å, which the peak height in the (gr) Ch + −FAs is less than g(r) Cl − −FAs (see Fig. 2b). It seems that the interaction between Cl − anion, and HBD leads to a decrease in the strong correlation between anion and cation. To investigate the structural properties of the binary mixtures with the different ratio of FAs and Ch + /Cl − salt, the RDFs was calculated at 353 K. For this purpose, the HA atom of FAs was chosen as a reference, and the distribution of the Cl − anions around them was investigated. A significant interaction is observed between HBA and HBD in the binary mixtures with the mole percent of CAP at 70.0% (see Fig. 2c). The coordination number (CN) was calculated to the number of anions around FAs in the binary mixtures with different ratios (see Table S1). The coordination number can be obtained from the first peak of the RDF curve between species 18 . For g(r) HBD_HBA , the coordination number is 0.1871 in the binary mixture with the mole percent of CAP at 30.0%, whereas, CN of RDF between HBA and HBD in the binary mixture with the mole percent of CAP at 70.0%, is 0.6432. Add the acid molecules to the simulation box, coordination numbers corresponding to the first solvation shell between HBA and HBD were increased in the binary mixture with the mole percent of CAP at 70.0%. In DES with 50% of CAP, the coordination number changes from 0.1871 to 0.3958. In contrast, DES with 50% of LUA has a coordination number of 0.2742. Decreasing CN between HBA and HBD can be attributed to the more favorable interaction of acid and salt in acid with the shorter chain length. It seems that caprylic acid molecules, which have a − 0.6269 charge on the OA atom, form a much stronger hydrogen bond with salt as compared to the LUA molecules, which have a − 0.6136 charge. RDFs between HBA and HBD were analyzed to characterize the structural correlation in aqueous solutions. In order to investigate the effect of water on the distribution of HBD around HBA, g(r) HBD−HBA was compared in pure and aqueous solutions of DES. Figure 2d shows RDFs between choline chloride salt and FAs molecules before and after mixing with water. It can be seen from Fig. 2d that the heights of peaks corresponding to g(r) HBA−HBD significantly decrease in water. The decreasing trend of the CN of g(r) between HBD and HBA is a confirmation of the decreasing distribution of HBD around HBA in the adjacent water (see Table S1). The coordination number related to RDF between HBA and HBD has also been used to determine structural changes in water. The structural properties of the binary mixtures may be greatly influenced by the adjacent water molecules (see Fig. 2d). The CN related to RDF between anion and HA atom of   Table S2 (Geometrical criteria were used to find H-bonds obtained in the RDF and the CDF). Hydrogen bonds were defined by a cut-off distance of less than 3.5 Å and a cut-off angle of less than180°. For obtaining the average number of H-bonds, Gaussian functions fitted to the number of H-bond data. www.nature.com/scientificreports/ where σ is the standard deviation, a and X represent the adjustable parameter and the average number of H-bond, respectively 19 . In this study, we focused on the hydrogen bonding interactions between the anion and the hydroxyl group of the cation and the -COOH group of FAs. The average number of H-bond between FAs molecules and salt is drawn in the binary mixtures with the different ratio of FAs and Ch + /Cl − salt (see Fig. 4a The relative percent of occupancy was analyzed to demonstrate the total time a unique hydrogen bond between two species 20 . The persistence of hydrogen bonds are obtains by using occupancy analysis 20 . Figure 5a shows the H atom of carboxylate group, HA, of FAs has a more stable interaction with anions of HBA. The stability of hydrogen bonds between acids and choline chloride salt were investigated at mixtures with different ratio of FAs (see Fig. 5b). Choline chloride salt in mixtures with 70% acid form strong and stable hydrogen bonds with the carboxylate of FAs. Hydrogen bond interaction is disturbed for all studied systems in the adjacent water, but the persistence of hydrogen bonds of the binary mixtures with 70% of FA is high as compared to other ratios (see Fig. 5c).

Nonbonded interaction energy and relative stability factor. Calculation of the non-bonded inter-
action energy is an important analysis for understanding intermolecular interactions in binary mixtures. The non-bonded energy includes terms such as short-range van der Waals (vdW) interactions and long-range electrostatic interactions which are computed Lennard-Jones (12-6) function and coulombic equation, respectively. The vdW and Coul interactions are calculated by the following equations: where q i andq j represent the atomic partial charges. D o and R o are the depth of the potential and the equilibrium atomic separation, respectively 21 .
The nonbonded interaction energies between different species of each system are listed in Table 2. In this section, a relative measure was computed to quantify the stability of DESs based on choline chloride and fatty acid in an aqueous solutions. To obtain the relative stability of hydrophobic DESs in water, the interaction energies between HBA-HBD, HBA-water, and HBD-water pairs were discussed in different ratios of fatty acids (caprylic acid and Lauric acid) and choline chloride. Relative stability is related to the stability of DES systems in water. Since the interaction energies between water and thymol are almost equal, it follows that the stability factor was a function of the HBD−water interaction . The relative stability factor (S) is defined as: To study the non-bonded interaction energy, Cl − −Water. H, Ch + −Water. H, and FAs. HA−Water. O interactions were selected due to the distributions of water molecules around carboxylate (-COO-) of FAs, hydroxyl (OH-) of Ch + cations, and chloride anion. E vdW and E Coul amounted to − 7.23 kcal/mol and − 116.044 kcal/ mol for HA. Fas-Cl − interaction in the binary mixture with 70% of FAs at the adjacent water, which E total is 31% more than the interaction energy of anions and CAP molecules in the binary mixture with 50% of FAs. To test the possibility that binary mixtures with a higher percentage of FAs composition are very stable in aqueous solutions, E total between HBA and HBD was compared for the binary mixtures with a molar percentage of 30%, 50 and 70% of FAs (ie. luaric acid and caprylic acid) at the adjacent of water. The results showed that the interaction between FAs molecules and Cl − anion was much more favorable in the binary mixture containing 70% of FAs than in other percentages.
It should be noted that, the value of the E total between LUA molecules and anion is than that of the LUA. H-Ch + .O interaction energy in the binary mixture at pure state. Furthermore, the E Coul contribution decreased an order of magnitude lower in LUA. H-Cl − interaction. In general, It seems that the choline chloride salt has a higher affinity toward the water-rich phase compared to fatty acids, so the interaction between the FAs molecules and the Cl − anion is dramatically reduced in the presence of water. However, the stability factor is a more suitable parameter to investigate the stability of binary mixtures in water. The stability factor values were calculated by www.nature.com/scientificreports/ using MD simulation (see Table 2). The stability factor values calculated on the basis of the simulated result is in range between 0.1 and 0.9. According to these results, the binary mixtures with a molar percentage of 70% of FAs have more stability in the adjacent water.
Spatial distribution functions. The spatial distribution function (SDF) of the species around each other was obtained through the TRAVIS package. The SDF analysis was used for understanding the microscopic structure of the binary mixtures in the adjacent water 22 . Green isosurfaces correspond with the spatial distribution of the anion around the FAs molecule, yellow isosurfaces are the distribution of the choline cation around reference www.nature.com/scientificreports/ FAs and red isosurfaces are the spatial distribution of fatty acid molecules around significant FAs in the binary mixtures with the mole percent at 50% of CAP (see Fig. 6(Panel a)). As can be appreciated from the colored SDFs www.nature.com/scientificreports/ for different species, the more active (-COOH) site of FAs molecules is surrounded by Cl − anions in the binary mixtures (see Fig. 6(Panel b)). The highly distributed anion around the FAs molecule can be seen in Fig. 6(Panel c). It should be noted that the anions are mainly located on the inner surface of the distribution of species around reference FAs molecules. Also, the distribution of acid molecules around choline is remarkable in the binary mixture.There are almost no significant changes in the distribution process of species with changing the number of acids from 300 to 500 in the box simulation. However, an isosurface shows that there is a rather high degree of association of HBA and HBD in the TLA with 70% LUA. Such a trend might appear reasonable since the spatial distribution of the anion around the cation was decreased in this ratio (see Fig. 6(Panel d)). SDFs reveal the fact that the distribution of species around each other is disturbed in water. The distribution of HBA around the central HBD gradually decreases as the binary mixture becomes more dilute (see Fig. 6(Panel c)). It seems that Cl − anions and Ch + cations prefer to be positioned near the H atoms of water molecules in the water/DES mixtures. However, by adding water molecules to the simulation box, the correlations between HBA and HBD were weak in the binary mixtures of 30% of CAP and 30% of LUA (see Fig. 6(Panel c)). The water molecules are evenly distributed around salt. However, It seems the distribution of water molecules around the FAs to be negligible in the binary mixtures with a molar percentage of 70% FAs (see Fig. 6(Panel d)).
Dynamic and transport properties. Rotational and lateral diffusion, mean square displacement. Accurate predictions of the dynamic behavior and the microscopic motion of the species will be interesting topics for future studies due to the widespread use of DESs in industry 23 . The dynamic properties of the eutectic solvents were also estimated from the self-diffusion coefficients of species by using the Einstein relation. Self-diffusion coefficients, D self , given by: where � − → r i (t + τ) − − → r i (t) 2 � stands for the mean-square displacement (MSD) of species i. The MSD is defined as The well-known β exponent was used to determine the location of the diffusive regime by Del Pópolo and Voth 24 . To determine the β parameter, the slope of MSD versus simulation time, was calculated as β = dlog 10 < �r(t) 2 > dlog 10 t www.nature.com/scientificreports/ According to the opinions of Del Pópolo and Voth 25 , D self was calculated from the MD simulations at the diffusive regime (β = 1) . Self-diffusion coefficients (D self ) of species were calculated from the slopes of MSD-time curves. The D self data are listed in Table 3. D self of Cl − anion is 0.2214, 0.4172 and 0.1927 A • 2 ns −1 or the binary www.nature.com/scientificreports/ mixtures with molar percentage of 30%, 50% and 70% CAP, respectively. D self of Cl − anion is 0.1927 A • 2 ns −1 in the binary mixtures with 70% of CAP, while the D self of Cl − corresponds to the same ratio in adjacent water is 124.2521 A • 2 ns −1 . When comparing the D self of species, we observed significant changes in the D self . The meansquare displacement (MSD) analysis is a very useful for tracking atomic motion. The mean square displacements (MSD) of species were plotted in the binary system versus the simulation time at 0−30 ns. The MSD functions at the binary mixture with a molar percentage of 70% LUA exhibit a hydrogen-bonding network that strongly limits the migration of molecules. There is one strong inermolecular O-H⋯Cl − hydrogen bond (FAs.H …Cl − ), resulting in a decrease in the slope of MSDs. These results are very good qualitative agreement with results of Combined Distribution Functions analysis (Fig. 7a,b).

Thermo-physical properties analysis.
To obtain the shear viscosity in molecular dynamics simulation, the shear autocorrelation function was calculated using the Green-Kubo method 26 . The shear viscosity is given by the following expression: where η is the shear viscosity, V and T represent the volume and the temperature of system, respectively. k B is the Boltzmann constant.P xy refers the off-diagonal element of the stress tensor 27 . After 50 ns NPT, five short runs, each of 1 ns in length, was performed with storing each frame. The shear viscosities for the studied systems are presented in Table 4. The shear viscosity of the binary mixture with 30% CAP at 353 K is 6.398 mPaS while adding FAs molecules increased the viscosity of the binary mixture at 70% CAP (see Table 4). Similarly, a small change in the shear viscosity was observed by changing the percentage of LUA from 30 to 50%. So, the binary mixtures with a molar percentage of 70% have the greatest shear viscosity. There is evidence of strong hydrogen bonding between Cl − anions of HBA and the headgroup of HBD at the binary mixtures with 70% FAs.
Velocity autocorrelation functions. The dynamics behavior of species in the binary mixtures can be studied by the normalized velocity autocorrelation function (VACF).C v can be calculated from the follows equation: where V c (t) related to the center of mass velocity of specie and the angular brackets is ensemble average over all time origins 28 .
A total of five simulation runs were performed for each of 1 ns in length. The input files of short runs were prepared from the latest step 50 ns in the ensemble NPT. To calculate VACFs, velocity data were collected in three directions in the x-axis, y-axis, and z-axis directions, with an interval of 10 ns. In order to confirm the statistical accuracy of the VACFs results and improve their reproducibility, an average of five runs was conducted for each of the studied systems. The computed VACFs of species for the binary mixture of FAs and [Ch + ] [Cl − ] with the molar percentage of 50% CAP are compared in Fig. S5. VACFs of species presented here show that Cl − anions velocities randomizes sooner than cation and FAs molecules. The point here is that the second zero is the velocity randomization time in the VACFs. The mean collision times (first zero) for Cl − anions at the molar percentage of 70% CAP are estimated at around 5 ps, whereas, the value of this quantity for a similar molar percentage is around 50 in the binary mixture of LUA and [Ch + ] [Cl − ] (Fig. 8a).The depth well of the first minimum of VACF for the binary mixtures with the molar percentage of 70% CAP confirms the high density in this ratio (see Table 4). The first depth well of VACF was decreased by adding water molecules to the simulation box. So that the negative region of the VACF cannot be clearly distinguished in water (Fig. 8b,c). The VRD(τ) is often defined as the normalized sum over the dot product between the vector at some time t, − → a i (t) , and the same vector at some later time t + τ, − → a i (t + τ ) for all starting times t 29 . To examine the reorientation of vectors, we have chosen the following bond vector. In the binary mixtures, the reorientation dynamics of the different vector of CAP molecules are visualized in Fig. 9a. The blue, red, and green lines show the reorientation of OA-HA bond vectors within each of the FAs molecules at the binary mixtures with the binary mixtures  www.nature.com/scientificreports/ with 30, 50 and 70% of CAP , respectively. It can be seen that the orientation of OA -HA bond vectors in the binary mixtures with 30 and 50% of FAs is faster than that of the binary mixtures with 70% of FAs (see Fig. 9b and Fig. S5). The first peak of the cation-FAs and anion-FAs RDFs in the binary mixtures with 70% of FAs is   Fig. 10a,b). It seems that the orientation of bond vectors is influenced by water molecules.

Conclusion
MD simulation was used to investigate the stability of DESs based on fatty acids and choline chloride in the adjacent water. For this purpose, the structural and dynamic properties of binary mixtures were investigated at different ratio of FAs and [Ch + ][Cl − ]. From the study on the structural properties of solvents, we infer that chloride anions play a key role in the formation of DESs. We mention that the structural properties were indicated   www.nature.com/scientificreports/ H-bond between the hydroxyl oxygen of cations and hydrogen atoms of carboxyl (-COOH) of FAs. Significant differences in intermolecular interaction of DESs were observed in the adjacent water. However, the binary mixtures containing 70% FAs, which anions have tend to form strong hydrogen bonds, are commonly the very stable mixture in water.